capture program drop cjive_simulations
program cjive_simulations
	syntax anything
	tokenize `anything'
	local B = `1' //number of iterations
	local N = `2' //number of observations
	local k = `3' //number of judges
	local c = `4' //number of clusters (max N, min k)
	local sigvc = `5'
	local rho = `6'
	local alpha = `7'
	local gamma = `8'

	************************************************
			*** Input Variables ***
	************************************************
	

	
	*Other Standard Values
	local sigvi = sqrt(1 - `sigvc'^2)
	local tau = 0.327
	
	*Matrix of Results
	mat results = J(`B',10,.) //Results Matrix
	mat colnames results = biasols seols bias2sls se2sls biasliml seliml biasjive sejive biasjkc sejkc
	tempfile results

	************************************************
			*** Data Creation ***
	************************************************

	*Different for each iteration
	forvalues iter = 1/`B'{
		clear
		scalar t1 = c(current_time)
		
		*Individual Level Data
		qui set obs `N'
		
		*Covariates
		gen x = invnormal(uniform())
		gen a = invnormal(uniform())
		
		*Judge Level Data
		gen p = `alpha' * x + a
		gen judgeid = ceil(`k'*normal(p / sqrt(1 + `alpha'^2)))

		*Deal with clusters
		gen cluster = ceil(`c' / `k' *judgeid + `c' / `k' *runiform()) - floor(`c' / `k')
		sort cluster
		gen vc = `sigvc' * invnormal(uniform())
		qui replace vc = vc[_n-1] if cluster[_n] == cluster[_n-1]
		gen vi = `sigvi' * invnormal(uniform())
		gen v = vc + vi

		*Everything else
		gen d = p + x + `gamma' * v //FIXME
		gen epsilon = `rho'*v+sqrt(1-`rho'^2)*invnormal(uniform())
		gen y = `tau' * d + x + epsilon
		
		

		************************************************
				*** Estimation ***
		************************************************
		
		*OLS
		qui reg y x d ,vce(cluster cluster)
		mat results[`iter',1]=_b[d]-`tau',_se[d]

		*2SLS
		qui ivregress 2sls y x (d = i.judge),vce(cluster cluster)
		mat results[`iter',3]=_b[d]-`tau',_se[d]

		*LIML
		qui ivregress liml y x (d = i.judge),vce(cluster cluster)
		mat results[`iter',5]=_b[d]-`tau',_se[d]

		*JIVE
		qui manyiv y x (d = i.judge),cluster(cluster)
		mat results[`iter',7]=e(b)[1,5]-`tau',e(se)[3,5]

		*CJIVE (Covariates Okay)
		qui xi i.judge
		qui cjive y x (d = _Ijudge*), cluster(cluster)
		mat results[`iter',9]=_b[d]-`tau',_se[d]

		scalar t2 = c(current_time)
		local iter_time = (clock(t2, "hms") - clock(t1, "hms")) / 1000
		if mod(`iter', 10) == 1 {
			disp ""
			disp "Completed `iter' out of `B', current iteration time: `iter_time' seconds"
			mat A = results[1..`iter' ,1..10]
			mat ones = J(`iter',1,1)
			mat sum = ones' * A
			mat means = sum / `iter'
			mat list means
		}

	}
	*Save Results
	clear
	svmat results,names(col)
	foreach type in ols 2sls liml jive jkc {
		gen covered`type' = bias`type'-1.96*se`type'<=0 & 0<=bias`type'+1.96*se`type' if bias`type'!=.
	}
	sum bias*
	sum se*
	sum covered*

	end
